.ipynb .pdf

repository open issue

Binder

Contents

Temperature, Salinity, and Stratification

import xarray as xr
from matplotlib import pyplot as plt
from IPython.display import SVG, display, Image, display_svg
import numpy as np
import gsw
%matplotlib inline
import holoviews as hv
from holoviews.operation.datashader import regrid
hv.extension('bokeh')
plt.rcParams['figure.figsize'] = (12,7)

# masks
woa_mask = xr.open_dataset('http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NODC/.WOA09/'
                             '.Masks/.basin/dods/').rename({'X':'lon', 'Y':'lat', 'Z': 'depth'})
basin_names = woa_mask.basin.attrs['CLIST'].split('\n')
basins_main = { n: basin_names[n-1] for n in [1,2,3,10] }
basins_marginal = { n: basin_names[n-1] for n in [4,5,6,7,11] }

# data
url_base = 'http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NODC/.WOA09/.Grid-1x1/.Annual'
url_temp = url_base + '/.temperature/.t_an/dods'
url_salt = url_base + '/.salinity/.s_an/dods '
woa_temp = xr.open_dataset(url_temp)
woa_temp_nobs = xr.open_dataarray(url_base + '/.temperature/.t_dd/dods')
woa_salt = xr.open_dataset(url_salt)
woa = xr.merge([woa_temp, woa_salt, woa_temp_nobs, woa_mask]).isel(time=0)
woa.load()

<xarray.Dataset>
Dimensions:  (depth: 33, lat: 180, lon: 360)
Coordinates:
  * lon      (lon) float32 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
  * depth    (depth) float32 0.0 10.0 20.0 30.0 ... 4000.0 4500.0 5000.0 5500.0
  * lat      (lat) float32 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
    time     datetime64[ns] 2008-01-01
Data variables:
    t_an     (depth, lat, lon) float32 nan nan nan nan nan ... nan nan nan nan
    s_an     (depth, lat, lon) float32 nan nan nan nan nan ... nan nan nan nan
    t_dd     (depth, lat, lon) float64 nan nan nan nan nan ... nan nan nan nan
    basin    (depth, lat, lon) float32 nan nan nan nan nan ... nan nan nan nan

pres = gsw.p_from_z(-woa.depth, 0)[:,None,None]
s_a = gsw.SA_from_SP(woa.s_an,
                     pres,
                     woa.lon.data[None,None,:],
                     woa.lat.data[None,:,None])
c_t = gsw.CT_from_t(s_a, woa.t_an, pres)
sigma_0 = gsw.sigma0(s_a, c_t)
sigma_2 = gsw.sigma2(s_a, c_t)
sigma_4 = gsw.sigma4(s_a, c_t)

for data, label in ((s_a, 's_a'), (c_t, 'c_t'),
                    (sigma_0, 'sigma_0'), (sigma_2, 'sigma_2'), (sigma_4, 'sigma_4')):
    woa[label] = xr.DataArray(data, dims=woa.t_an.dims, coords=woa.t_an.coords)

woa

<xarray.Dataset>
Dimensions:  (depth: 33, lat: 180, lon: 360)
Coordinates:
  * lon      (lon) float32 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
  * depth    (depth) float32 0.0 10.0 20.0 30.0 ... 4000.0 4500.0 5000.0 5500.0
  * lat      (lat) float32 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
    time     datetime64[ns] 2008-01-01
Data variables:
    t_an     (depth, lat, lon) float32 nan nan nan nan nan ... nan nan nan nan
    s_an     (depth, lat, lon) float32 nan nan nan nan nan ... nan nan nan nan
    t_dd     (depth, lat, lon) float64 nan nan nan nan nan ... nan nan nan nan
    basin    (depth, lat, lon) float32 nan nan nan nan nan ... nan nan nan nan
    s_a      (depth, lat, lon) float64 nan nan nan nan nan ... nan nan nan nan
    c_t      (depth, lat, lon) float64 nan nan nan nan nan ... nan nan nan nan
    sigma_0  (depth, lat, lon) float64 nan nan nan nan nan ... nan nan nan nan
    sigma_2  (depth, lat, lon) float64 nan nan nan nan nan ... nan nan nan nan
    sigma_4  (depth, lat, lon) float64 nan nan nan nan nan ... nan nan nan nan

%%opts Image [logz=True width=700 height=400 colorbar=True toolbar='above']
hv.Dataset(woa.t_dd).to(hv.Image, ['lon', 'lat'], label='Number of Observations')

%%opts Image [width=700 height=400 colorbar=True toolbar='above' tools=['hover']] (cmap='Category20')
hv.Dataset(woa.basin.where(woa.basin<20)).to(hv.Image, ['lon', 'lat'], label='Basin')

def map_section_plot(data, label, data_range, levels, cmap='viridis'):
    ds = hv.Dataset(data)
    img = ds.to(hv.Image, ['lon', 'lat'], label=label)
    img = img.redim.range(**{data.name: data_range})

    posx = hv.streams.PointerX()
    vline = hv.DynamicMap(lambda x: hv.VLine(x or 180.5), streams=[posx])

    lats = np.arange(-90, 91)
    def cross_section(x):
        mesh = hv.QuadMesh((lats, woa.depth, data.sel(lon=x if x else 180, method='nearest').data[:-1]),
                           kdims=['lat', 'depth'], label=(label + ' Section'))
        contours = hv.operation.contours(mesh, levels=levels,
                                         overlaid=False, filled=True)
        return contours

    
    crosssection = hv.DynamicMap(cross_section, streams=[posx])
    layout = (img * vline) + crosssection 
    
    # set view options
    dict_spec = {'Image': {'plot': dict(colorbar=True, toolbar='above', width=600),
                           'style': dict(cmap=cmap)},
                 'NdOverlay': {'plot': dict(bgcolor='grey', invert_yaxis=True, width=600, height=200)},
                 'Layout': {'plot': dict(show_title=False)},
                 'Polygons': {'plot': dict(colorbar=True, invert_yaxis=True, width=600),
                              'style': dict(cmap=cmap, line_width=0.1)}
                 }
    layout.cols(1)
    layout = layout.opts(dict_spec)   
    return layout

map_section_plot(woa.t_an, 'Temperature', (-2,30), np.arange(-2,31), cmap='RdBu_r')

map_section_plot(woa.s_an, 'Salinity', (32,38), np.arange(32,38.1, 0.25))

map_section_plot(woa.sigma_0, 'Potential Density (sigma0)', (22.,28.), np.arange(22,28.1, 0.25))

woa['thick'] = woa.depth.diff(dim='depth')
woa['thick'][0] = 5
# add an area element
earth_radius =  6.371e6
tot_area = (4 * np.pi * earth_radius**2)

woa['area'] = (earth_radius**2 *
               np.radians(1.) * np.radians(1.) *
               np.cos(np.radians(woa.lat)))
t, s, dz, a, basin = xr.broadcast_arrays(
                *xr.align(
                    woa['t_an'], woa['s_an'],
                    woa['thick'], woa['area'], woa['basin']
                )
              )
vol = dz*a

/Users/rpa/miniconda3/envs/geo_scipy/lib/python3.6/site-packages/ipykernel_launcher.py:13: DeprecationWarning: xarray.broadcast_arrays is deprecated: use xarray.broadcast instead
  del sys.path[0]

woa['vol'] = woa.thick * woa.area

fig, (ax1, ax2) = plt.subplots(2)

ax1.hist(t.where(t).values.ravel(), bins=100, normed=True, range=(-2,30),
         weights=vol.where(t).values.ravel());
ax2.hist(s.where(s).values.ravel(), bins=100, normed=True, range=(32,38),
         weights=vol.where(s).values.ravel());
ax1.set_xlabel(r'Temperature ($^\circ$C)')
ax2.set_xlabel("Salinity (g / kg)")
plt.tight_layout()

/Users/rpa/miniconda3/envs/geo_scipy/lib/python3.6/site-packages/numpy/lib/function_base.py:748: RuntimeWarning: invalid value encountered in greater_equal
  keep = (tmp_a >= mn)
/Users/rpa/miniconda3/envs/geo_scipy/lib/python3.6/site-packages/numpy/lib/function_base.py:749: RuntimeWarning: invalid value encountered in less_equal
  keep &= (tmp_a <= mx)
/Users/rpa/miniconda3/envs/geo_scipy/lib/python3.6/site-packages/numpy/lib/function_base.py:748: RuntimeWarning: invalid value encountered in greater_equal
  keep = (tmp_a >= mn)
/Users/rpa/miniconda3/envs/geo_scipy/lib/python3.6/site-packages/numpy/lib/function_base.py:749: RuntimeWarning: invalid value encountered in less_equal
  keep &= (tmp_a <= mx)

Temperature and Salinity Distribution

fig

from matplotlib.colors import LogNorm
fig, ax = plt.subplots()
hb = ax.hexbin(s.values.ravel(), t.values.ravel(),
           C=vol.values.ravel(), reduce_C_function=np.sum,
           extent=(30,40,-3,30), gridsize=50, bins='log',
           cmap='RdYlBu_r')
plt.colorbar(hb)
ax.set_ylabel(r'Temperature ($^\circ$C)')
ax.set_xlabel("Salinity (g / kg)")
hb.set_clim([12,18])

Joint T/S Distribution

fig

#pal = seaborn.palettes.color_palette(n_colors=len(basins_main))
fig, ax = plt.subplots()
for (bnum, bname) in basins_main.items():
    ax.plot(s.where(basin==bnum).values.ravel(), t.where(basin==bnum).values.ravel(),
             '.',  alpha=0.02, label=bname)
ax.set_xlim(30,38)
ax.legend(loc='upper left', frameon=True)
# blue green red purple

<matplotlib.legend.Legend at 0x10385a3c8>

T/S Distribution by Basin

fig

Image("02_images/THC_Ballarotta2014.png", width=600)

# the goods are here
# https://www.nodc.noaa.gov/OC5/3M_HEAT_CONTENT/
# http://data.nodc.noaa.gov/thredds/catalog/woa/heat_content/heat_content/catalog.html
url = ('https://data.nodc.noaa.gov/thredds/dodsC/woa/heat_content/'
       'heat_content/heat_content_anomaly_0-700_pentad.nc')
hc = xr.open_dataset(url, decode_times=False).load()
hc

<xarray.Dataset>
Dimensions:             (depth: 1, lat: 180, lon: 360, nbounds: 2, time: 59)
Coordinates:
  * lat                 (lat) float32 -89.5 -88.5 -87.5 -86.5 ... 87.5 88.5 89.5
  * lon                 (lon) float32 -179.5 -178.5 -177.5 ... 177.5 178.5 179.5
  * time                (time) float32 30.0 42.0 54.0 66.0 ... 702.0 714.0 726.0
Dimensions without coordinates: depth, nbounds
Data variables:
    crs                 int32 -2147483647
    lat_bnds            (lat, nbounds) float32 -90.0 -89.0 -89.0 ... 89.0 90.0
    lon_bnds            (lon, nbounds) float32 -180.0 -179.0 ... 179.0 180.0
    depth_bnds          (depth, nbounds) float32 0.0 700.0
    climatology_bounds  (time, nbounds) float32 0.0 60.0 12.0 ... 696.0 756.0
    h18_hc              (time, depth, lat, lon) float32 nan nan ... 0.036448
    pent_h22_WO         (time) float32 -6.28869 -4.89276 ... 13.3964 14.1427
    pent_h22_se_WO      (time) float32 1.84269 1.83688 ... 0.134667 0.131762
    pent_h22_NH         (time) float32 -1.43756 -0.659048 ... 6.14617 6.52268
    pent_h22_se_NH      (time) float32 0.760144 0.708205 ... 0.061274 0.059592
    pent_h22_SH         (time) float32 -4.85115 -4.2337 ... 7.25029 7.62019
    pent_h22_se_SH      (time) float32 1.08254 1.12867 ... 0.073392 0.072169
    pent_h22_AO         (time) float32 -2.60102 -2.40377 ... 5.47195 5.70196
    pent_h22_se_AO      (time) float32 0.549281 0.451574 ... 0.042949 0.042403
    pent_h22_NA         (time) float32 -0.971802 -0.658995 ... 3.23728 3.40875
    pent_h22_se_NA      (time) float32 0.323197 0.22556 ... 0.024828 0.024125
    pent_h22_SA         (time) float32 -1.62921 -1.74478 ... 2.23473 2.29322
    pent_h22_se_SA      (time) float32 0.226085 0.226015 ... 0.018121 0.018279
    pent_h22_PO         (time) float32 -2.50564 -1.25884 ... 4.29607 4.88475
    pent_h22_se_PO      (time) float32 0.852305 0.907993 ... 0.063972 0.061715
    pent_h22_NP         (time) float32 -0.396412 0.065448 ... 2.13457 2.36667
    pent_h22_se_NP      (time) float32 0.363998 0.411451 ... 0.031127 0.03018
    pent_h22_SP         (time) float32 -2.10923 -1.32429 ... 2.16149 2.51809
    pent_h22_se_SP      (time) float32 0.488307 0.496542 ... 0.032845 0.031536
    pent_h22_IO         (time) float32 -1.19054 -1.23783 ... 3.61581 3.54261
    pent_h22_se_IO      (time) float32 0.440468 0.476714 ... 0.026263 0.026152
    pent_h22_NI         (time) float32 -0.077804 -0.073198 ... 0.76171 0.733754
    pent_h22_se_NI      (time) float32 0.072317 0.070599 ... 0.003836 0.003797
    pent_h22_SI         (time) float32 -1.11274 -1.16463 ... 2.85409 2.80886
    pent_h22_se_SI      (time) float32 0.368151 0.406115 ... 0.022427 0.022355
    basin_mask          (lat, lon) float64 nan nan nan nan ... 1.0 1.0 1.0 1.0
Attributes:
    Conventions:                     CF-1.6
    title:                           Ocean Heat Content anomalies from WOA09 ...
    summary:                         Mean ocean variable anomaly from in situ...
    references:                      Levitus, S., J. I. Antonov, T. P. Boyer,...
    institution:                     National Oceanographic Data Center(NODC)
    comment:                         
    id:                              heat_content_anomaly_0-700_pentad.nc
    naming_authority:                gov.noaa.nodc
    time_coverage_start:             1955-01-01
    time_coverage_duration:          P63Y
    time_coverage_resolution:        P05Y
    geospatial_lat_min:              -90.0
    geospatial_lat_max:              90.0
    geospatial_lon_min:              -180.0
    geospatial_lon_max:              180.0
    geospatial_vertical_min:         0.0
    geospatial_vertical_max:         700.0
    geospatial_lat_units:            degrees_north
    geospatial_lat_resolution:       1.00 degrees
    geospatial_lon_units:            degrees_east
    geospatial_lon_resolution:       1.00 degrees
    geospatial_vertical_units:       m
    geospatial_vertical_resolution:  
    geospatial_vertical_positive:    down
    creator_name:                    Ocean Climate Laboratory
    creator_email:                   NODC.Services@noaa.gov
    creator_url:                     http://www.nodc.noaa.gov
    project:                         World Ocean Database
    processing_level:                processed
    keywords:                        <ISO_TOPIC_Category> Oceans</ISO_TOPIC_C...
    keywords_vocabulary:             ISO 19115
    standard_name_vocabulary:        CF-1.6
    contributor_name:                Ocean Climate Laboratory
    contributor_role:                Calculation of anomalies
    featureType:                     Grid
    cdm_data_type:                   Grid
    nodc_template_version:           NODC_NetCDF_Grid_Template_v1.0
    date_created:                    2018-01-18 
    date_modified:                   2018-01-18 
    publisher_name:                  US NATIONAL OCEANOGRAPHIC DATA CENTER
    publisher_url:                   http://www.nodc.noaa.gov/
    publisher_email:                 NODC.Services@noaa.gov
    license:                         These data are openly available to the p...
    Metadata_Conventions:            Unidata Dataset Discovery v1.0
    metadata_link:                   http://www.nodc.noaa.gov/OC5/3M_HEAT_CON...

# fix time to something xray likes
hc['time'] *= (365/12.)
tunits = hc.time.attrs['units'].replace('months', 'days')
hc['time'].attrs['units'] = tunits
hc = xr.decode_cf(hc)

fig, ax = plt.subplots()
hc['h18_hc'].mean(dim=('lon','lat','depth')).plot(ax=ax)
ax.set_ylabel(hc.h18_hc.attrs['units'])
plt.title('Ocean Heat Content Anomaly')
plt.grid()

fig

import datashader as ds
from holoviews.operation.datashader import datashade, dynspread

df = woa[['s_a', 'c_t', 'vol', 'basin']].to_dataframe().reset_index()
df = df.dropna()
df = df[df['basin'] <= 20]
df['basin'] = df['basin'].astype('category')
hv_df = hv.Dataset(df)
hv_df

:Dataset   [depth,lat,lon,s_a,c_t,vol,basin,time]

df.basin.unique()

[10.0, 1.0, 3.0, 2.0, 7.0, ..., 16.0, 17.0, 18.0, 19.0, 20.0]
Length: 20
Categories (20, float64): [10.0, 1.0, 3.0, 2.0, ..., 17.0, 18.0, 19.0, 20.0]

points = hv.Points(hv_df, kdims=['s_a', 'c_t'], vdims=['basin', 'vol'])
print(points)

:Points   [s_a,c_t]   (basin,vol)

%%opts RGB [width=400 height=300 colorbar=True tools=['hover']]
datashade(points) + datashade(points, aggregator=ds.sum('vol'))

%%opts RGB [width=400 height=300 colorbar=True tools=['hover']]
dynspread(datashade(points, aggregator=ds.sum('vol'))) + datashade(points, aggregator=ds.sum('vol'))

%%opts RGB [width=600 height=400 colorbar=True]
datashade(points, aggregator=ds.count_cat('basin'))

previous

The Water Column

next

Air Sea Exchange